Kernel Author:
Bhishan Poudel, Data Scientist, Ph.D Astrophysics .
Bhishan Poudel, Data Scientist, Ph.D Astrophysics .
Data Description¶
Reference: https://www.kaggle.com/c/web-traffic-time-series-forecasting/data
Original data: train_1.csv
-----------------------------
rows = 145,063
columns = 551
first column = Page
date columns = 2015-07-01, 2015-07-02, ..., 2016-12-31 (550 columns)
file size: 284.6 MB
Data for modelling:
--------------------------------------------------------------------
Timeseries : Now You See Me es (Spain, random_state=42)
For ARIMA : we have only one timeseries (one column)
For sklearn : For linear regressor, ensember learners we can have many columns
For fbprophet: we have only dataframe with columns ds and y (additional cap and floor)Prophet Description¶
References:
- https://www.analyticsvidhya.com/blog/2018/05/generate-accurate-forecasts-facebook-prophet-python-r/
- https://www.kaggle.com/rpbenatti/dh-desafio-4-benatti
- https://www.kaggle.com/conorm97/prediction-using-prophet-with-us-holidays
- https://www.kaggle.com/springrid/prediction-prophet-movingavg-and-apple-vs-facebook
- https://www.kaggle.com/cpmpml/smape-weirdness
We use a decomposable time series model with three main model components: trend, seasonality, and holidays. They are combined in the following equation:
$$ y(t)=g(t)+s(t)+h(t)+\epsilon_{t} $$- g(t): piecewise linear or logistic growth curve for modelling non-periodic changes in time series
- s(t): periodic changes (e.g. weekly/yearly seasonality)
- h(t): effects of holidays (user provided) with irregular schedules
- εt: error term accounts for any unusual changes not accommodated by the model
Using time as a regressor, Prophet is trying to fit several linear and non linear functions of time as components.
Modeling seasonality as an additive component is the same approach taken by exponential smoothing in Holt-Winters technique .
We are, in effect, framing the forecasting problem as a curve-fitting exercise rather than looking explicitly at the time based dependence of each observation within a time series.
Trend parameters:
| Parameter | Description |
|---|---|
| growth | linear’ or ‘logistic’ to specify a linear or logistic trend |
| changepoints | List of dates at which to include potential changepoints (automatic if not specified) |
| n_changepoints | If changepoints in not supplied, you may provide the number of changepoints to be automatically included |
| changepoint_prior_scale | Parameter for changing flexibility of automatic changepoint selection |
Seasonality & Holiday Parameters:
| Parameter | Description |
|---|---|
| yearly_seasonality | Fit yearly seasonality |
| weekly_seasonality | Fit weekly seasonality |
| daily_seasonality | Fit daily seasonality |
| holidays | Feed dataframe containing holiday name and date |
| seasonality_prior_scale | Parameter for changing strength of seasonality model |
| holiday_prior_scale | Parameter for changing strength of holiday model |
In [1]:
# from fbprophet import Prophet
# help(Prophet)
Imports¶
In [1]:
import numpy as np
import pandas as pd
import seaborn as sns
sns.set(color_codes=True)
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
sns.set(context='notebook', style='whitegrid', rc={'figure.figsize': (12,8)})
plt.style.use('ggplot') # better than sns styles.
matplotlib.rcParams['figure.figsize'] = 12,8
import os
import time
# random state
random_state=100
np.random.seed(random_state)
# Jupyter notebook settings for pandas
#pd.set_option('display.float_format', '{:,.2g}'.format) # numbers sep by comma
from pandas.api.types import CategoricalDtype
np.set_printoptions(precision=3)
pd.set_option('display.max_columns', 100)
pd.set_option('display.max_rows', 100) # None for all the rows
pd.set_option('display.max_colwidth', 200)
import IPython
from IPython.display import display, HTML, Image, Markdown
print([(x.__name__,x.__version__) for x in [np, pd,sns,matplotlib]])
In [2]:
%%javascript
IPython.OutputArea.auto_scroll_threshold = 9999;
In [3]:
from sklearn.preprocessing import StandardScaler
In [4]:
from fbprophet import Prophet
pd.plotting.register_matplotlib_converters() # prophet needs this
from fbprophet.plot import add_changepoints_to_plot
In [5]:
from datetime import date
import holidays
In [6]:
from fbprophet.plot import plot_plotly
import plotly.offline as py
py.init_notebook_mode()
In [7]:
from fbprophet.diagnostics import cross_validation
In [8]:
import plotly
import plotly.offline as py
import plotly.graph_objs as go
import plotly.figure_factory as ff
import plotly.tools as tls
from plotly.tools import make_subplots
from plotly.offline import plot, iplot, init_notebook_mode
init_notebook_mode(connected=False)
[(x.__name__,x.__version__) for x in [plotly]]
Out[8]:
Useful Scripts¶
In [9]:
%%writefile ../models/util_prophet.py
import numpy as np
import pandas as pd
import seaborn as sns
sns.set(color_codes=True)
import matplotlib
import matplotlib.pyplot as plt
plt.style.use('ggplot') # better than sns styles.
matplotlib.rcParams['figure.figsize'] = 12,8
import plotly
import plotly.offline as py
import plotly.graph_objs as go
import plotly.figure_factory as ff
import plotly.tools as tls
from plotly.tools import make_subplots
from plotly.offline import plot, iplot, init_notebook_mode
from numba import jit
import math
# https://www.kaggle.com/cpmpml/smape-weirdness
@jit
def smape(y_true, y_pred):
A = y_true.to_numpy().ravel()
F = y_pred.to_numpy().ravel()[:len(A)]
return ( 200.0/len(A) * np.sum( np.abs(F - A) /
(np.abs(A) + np.abs(F) + np.finfo(float).eps))
)
# https://www.kaggle.com/c/web-traffic-time-series-forecasting/discussion/37232
@jit
def smape_fast(y_true, y_pred):
"""Fast implementation of SMAPE.
Parameters
-------------
y_true: numpy array with no NaNs and non-negative
y_pred: numpy array with no NaNs and non-negative
Returns
-------
out : float
"""
out = 0
for i in range(y_true.shape[0]):
if (y_true[i] != None and np.isnan(y_true[i]) == False):
a = y_true[i]
b = y_pred[i]
c = a+b
if c == 0:
continue
out += math.fabs(a - b) / c
out *= (200.0 / y_true.shape[0])
return out
def safe_median(s):
return np.median([x for x in s if ~np.isnan(x)])
def mean_absolute_percentage_error(y_true, y_pred):
"""Calculate the MAPE.
"""
return np.mean(np.abs((y_true - y_pred) / y_true)) * 100
def plot_actual_forecast_mpl(df, forecast_str_lst, forecast_lst):
"""Plot prophet forcast.
Parameters
-----------
df -- dataframe with columns ds,y (cap and floor are optional)
forecast_str_lst -- list of strings
forecast_lst -- list of forecasts
Example
--------
forecast_str_lst = ['forecast1', 'forecast2', 'forecast3','forecast4']
forecast_lst = [eval(i) for i in forecast_str_lst]
plot_actual_forecast_mpl(df, ['forecast1'], [forecast1])
"""
import matplotlib.pyplot as plt
plt.style.use('ggplot')
plt.figure(figsize=(12,8))
plt.plot(df.y,'b',label='original')
colors10_hex = ['#b03060','#ff0000', '#ff00ff',
'#67ceab', '#63c56c', '#225e31',
'#29b6f6', '#6495ed','#00008b',
'#ffa500']
for i,f in enumerate(forecast_str_lst):
forecast = forecast_lst[i]
plt.plot(forecast.yhat,c=colors10_hex[i],label=f)
plt.legend()
def plot_actual_forecast_sns(df, forecast_str_lst, forecast_lst):
"""Plot prophet forcast.
Parameters
-----------
df -- dataframe with columns ds,y (cap and floor are optional)
forecast_str_lst -- list of strings
Example
--------
forecast_str_lst = ['forecast1', 'forecast2', 'forecast3','forecast4']
forecast_lst = [eval(i) for i in forecast_str_lst]
plot_actual_forecast_sns(df2, forecast_str_lst, forecast_lst)
"""
import seaborn as sns
import matplotlib.pyplot as plt
plt.style.use('ggplot')
plt.figure(figsize=(12,8))
df_plot = df2[['y']]
df_plot.index = pd.to_datetime(df2['ds'])
for i,f in enumerate(forecast_str_lst):
forecast = forecast_lst[i]
ts = forecast.yhat
ts.index = pd.to_datetime(forecast.ds)
df_tmp = pd.DataFrame({f: ts})
df_plot = pd.concat([df_plot,df_tmp], axis=1)
sns.lineplot(data=df_plot)
def plot_actual_forecast_plotly(df, forecast_str_lst, forecast_lst):
"""Plot prophet forcast.
Parameters
-----------
df -- dataframe with columns ds,y (cap and floor are optional)
forecast_str_lst -- list of strings
forecast_lst -- list of forecasts
Example
--------
forecast_str_lst = ['forecast1', 'forecast2', 'forecast3','forecast4']
forecast_lst = [eval(i) for i in forecast_str_lst]
plot_actual_forecast_plotly(df2, forecast_str_lst,forecast_lst)
"""
from plotly.offline import plot, iplot, init_notebook_mode
df_plot = df[['y']]
df_plot.index = pd.to_datetime(df['ds'])
for i,f in enumerate(forecast_str_lst):
forecast = forecast_lst[i]
ts = forecast.yhat
ts.index = pd.to_datetime(forecast.ds)
df_tmp = pd.DataFrame({f: ts})
df_plot = pd.concat([df_plot,df_tmp], axis=1)
iplot([{'x': df_plot.index,'y': df_plot[col],'name': col}
for col in df_plot.columns
])
def plot_deltas(m):
"""Plot model params delta as bar plots.
Notes:
-------
1. If the barplot is all incresing downward,
we may need to change these quantities:
- changepoint_range=0.8 (default is 0.8)
- changepoint_prior_scale=0.7 (default is 0.05)
"""
import matplotlib.pyplot as plt
plt.style.use('ggplot')
deltas = m.params['delta'].mean(axis=0)
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111)
ax.bar(range(len(deltas)), deltas)
ax.set_ylabel('Rate change (delta)')
ax.set_xlabel('Potential changepoint')
fig.tight_layout()
def plot_deltas_plotly(m):
"""Plot prophet forecast params delta values.
"""
import plotly.graph_objs as go
from plotly.offline import plot, iplot
# data to plot
x = list(range(len(m.params['delta'][0])))
y = m.params['delta'].ravel().tolist()
# trace
trace = go.Bar(x= x,y=y,name='Change Points')
data = [trace]
fig = go.Figure(data=data)
iplot(fig)
def outliers_to_na(ts, devs):
"""Replace the outliers by na.
Then we can again fill na by 0.
Here, in this wikipedia data nans are given 0.
"""
median= ts['y'].median()
std = np.std(ts['y'])
for x in range(len(ts)):
val = ts['y'][x]
if (val < median - devs * std or val > median + devs * std):
ts['y'][x] = None
return ts
def convert_ts_to_prophet_df(ts):
"""Convert timeseries to dataframe required by prophet.
Parameters
-----------
ts: timeseries with index as datetime and have values
"""
df = pd.DataFrame(columns=['ds','y'])
df['ds'] = pd.to_datetime(ts.index)
df.index = df['ds']
df['y'] = ts.to_numpy()
return df
def remove_negs_from_forecast(forecast):
"""Replace negative forecasts by 0.
Parameters
----------
forecast -- dataframes returned by prophet
Example
--------
m1 = Prophet()
m1.fit(df1);
future1 = m1.make_future_dataframe(periods=60)
forecast1 = m1.predict(future1)
"""
forecast = forecast.copy()
forecast['yhat'] = forecast['yhat'].clip_lower(0)
forecast['yhat_lower'] = forecast['yhat_lower'].clip_lower(0)
forecast['yhat_upper'] = forecast['yhat_upper'].clip_lower(0)
return forecast
In [10]:
import sys
sys.path.append('../models')
from util_prophet import smape, smape_fast
from util_prophet import plot_actual_forecast_plotly, plot_deltas_plotly
Load the data¶
In [11]:
%%time
df_raw = pd.read_csv('../../data/wiki/train_1.csv',encoding='latin-1')
print(df_raw.shape) # (145063, 551) we have 145k data and 551 columns
df_raw.head()
Data Preprocessing¶
Data Sample¶
In [12]:
df_sample = df_raw.sample(5, random_state=42) # I like 100 but use 42 here.
df_sample
Out[12]:
In [13]:
ts = df_sample.iloc[1,1:].astype(np.float64)
ts.name = df_sample.iloc[1,0]
ts.index = pd.to_datetime(ts.index,format='%Y-%m-%d')
ts.head()
Out[13]:
In [14]:
pd.plotting.register_matplotlib_converters() # prophet needs this
ts.plot()
Out[14]:
In [15]:
iplot([{'x': ts.index,'y': ts.to_numpy()} ])
Modelling: prophet¶
Create dataframe with two columns: ds and y¶
In [16]:
# prophet expects two columns: ds and y
df1 = pd.DataFrame({'ds': ts.index, 'y': ts.to_numpy()})
print(df1.dtypes)
df1.head()
Out[16]:
In [17]:
# df1.plot(kind='line',x='ds',y='y')
In [18]:
iplot([{'x': df1['ds'], 'y': df1['y']}])
model1: default parameters¶
In [19]:
m1 = Prophet()
m1.fit(df1);
In [23]:
# future
future1 = m1.make_future_dataframe(periods=60)
print(df1.shape)
print(future1.shape)
future1.head().append(future1.tail())
Out[23]:
In [21]:
# forecast
forecast1 = m1.predict(future1)
print(forecast1.shape)
forecast1.head(2)
Out[21]:
In [23]:
forecast1.columns
Out[23]:
In [24]:
forecast1[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].tail()
Out[24]:
In [25]:
fig1 = m1.plot(forecast1)
fig2 = m1.plot_components(forecast1)
In [26]:
plot_actual_forecast_plotly(df1, ['forecast1'],[forecast1])
In [27]:
forecast1.tail()['yhat']
Out[27]:
In [28]:
# observation: this prediction only captures the trend
# problem : the trend is going down and predicts negative values
#
# todo:
# use logistic growth and give cap and floor for dataframes df and forecast.
Model2: saturation cap and floor¶
In [29]:
# set the cap to values of trend
df2 = df1.copy()
df2['cap'] = 500
df2['floor'] = 0.0
future2 = future1.copy()
future2['cap'] = 500
future2['floor'] = 0.0
# take growth = 'logistic'
m2 = Prophet(growth='logistic')
forecast2 = m2.fit(df2).predict(future2)
# Plotting both the forecast predictions and components
fig1 = m2.plot(forecast2)
fig2 = m2.plot_components(forecast2)
In [30]:
# Conclusion: in this case the prediction trend has max 500 and min 0.0
# as we have set the values.
Model3: Seasonality¶
In [31]:
# keep default growth='linear' instead of logistic
m3 = Prophet(growth='linear',
daily_seasonality=True,
weekly_seasonality=True,
yearly_seasonality=True)
forecast3 = m3.fit(df2).predict(future2)
fig1 = m3.plot(forecast3)
fig2 = m3.plot_components(forecast3)
In [32]:
plot_actual_forecast_plotly(df2, ['forecast3'],[forecast3])
In [33]:
forecast3.tail()['yhat']
Out[33]:
In [34]:
# observations:
# 1. trend goes downwards in the end and gives negative values
# 2. gives better fit than last time
# 3. peaks are captured better.
Model4: Changepoints¶
In [35]:
from fbprophet.plot import add_changepoints_to_plot
fig = m3.plot(forecast3)
a = add_changepoints_to_plot(fig.gca(), m3, forecast3)
In [36]:
# observation:
# 1. default changepoints = 0.8 takes only 80% of data to find changepoints
# 2. we see some peaks after 80% of data, increase changepoints
In [37]:
m4 = Prophet(daily_seasonality=True,
weekly_seasonality=True,
yearly_seasonality=True,
changepoint_range=0.9)
forecast4 = m4.fit(df2).predict(future2)
fig = m4.plot(forecast4)
a = add_changepoints_to_plot(fig.gca(), m4, forecast4)
In [38]:
plot_deltas_plotly(m4)
In [39]:
# observation:
# 1. when changepoint range increase, predictions becomes more negative
# upto about 17, then becomes more positive but still negative.
# 2. I can fiddle with default changepoint 0.8.
#
# todo:
# default changepoint prior scale = 0.05, increase it to make more flexible trend.
In [40]:
m5 = Prophet(daily_seasonality=True,
weekly_seasonality=True,
yearly_seasonality=True,
changepoint_range=0.8,
changepoint_prior_scale=0.7)
forecast5 = m5.fit(df2).predict(future2)
fig = m5.plot(forecast5)
a = add_changepoints_to_plot(fig.gca(), m5, forecast5)
In [41]:
plot_deltas_plotly(m5)
In [42]:
forecast_str_lst = ['forecast1', 'forecast2', 'forecast3','forecast4','forecast5']
forecast_lst = [eval(i) for i in forecast_str_lst]
plot_actual_forecast_plotly(df2, forecast_str_lst,forecast_lst)
In [43]:
# observations:
# 1. forecast5 is almost similar to 4 and 3
# 2. forecast5 gives better peaks than 3
Model6: Holidays¶
In [44]:
df_sample
Out[44]:
In [45]:
# observation:
# 1. I have chosen the wikipedia page of Now you see me in Spain as our
# time series dataframe. I will look at holidays in Spain.
In [46]:
from datetime import date
import holidays
In [47]:
# Select country
es_holidays = holidays.Spain(years = [2015,2016,2017])
es_holidays = pd.DataFrame.from_dict(es_holidays, orient='index')
es_holidays = pd.DataFrame({'holiday': 'Spain', 'ds': es_holidays.index})
print(es_holidays.shape)
es_holidays.head()
Out[47]:
In [48]:
m6 = Prophet(growth='linear',
daily_seasonality=True,
weekly_seasonality=True,
yearly_seasonality=True,
changepoint_range=0.8,
changepoint_prior_scale=0.7,
holidays=es_holidays)
# add holidays
m6.add_country_holidays(country_name='ES')
# forecast
forecast6 = m6.fit(df2).predict(future2)
# plot forecasts
fig1 = m6.plot(forecast6)
fig2 = m6.plot_components(forecast6)
Plotly Visualizations for prophet¶
In [49]:
from fbprophet.plot import plot_plotly
import plotly.offline as py
py.init_notebook_mode()
fig = plot_plotly(m6, forecast6)
py.iplot(fig)
SMAPE predictions¶
The formula for SMAPE (Symmetric Mean Absolute Percentage Error) is given below:
$$ S M A P E=\frac{100 \%}{n} \sum_{t=1}^{n} \frac{\left|F_{t}-A_{t}\right|}{\left(\left|A_{t}\right|+\left|F_{t}\right|\right) / 2} $$Where, F is forecast and A is the actual value of time series at given time t.
Python implementation:
def smape(A, F):
F = A[:len(A)]
return ( 200.0/len(A) * np.sum( np.abs(F - A) /
(np.abs(A) + np.abs(F) + np.finfo(float).eps))
)
Despite the name Symmetric, the smape is not actually symmetric. Take this example from wikipedia for an example:
The SMAPE is not symmetric since over- and under-forecasts are not treated equally. This is illustrated by the following example by applying the SMAPE formula:
Over-forecasting: At = 100 and Ft = 110 give SMAPE = 4.76%
Under-forecasting: At = 100 and Ft = 90 give SMAPE = 5.26%.In [50]:
def smape_np(A, F):
F = F[:len(A)]
return ( 200.0/len(A) * np.sum( np.abs(F - A) /
(np.abs(A) + np.abs(F) + np.finfo(float).eps))
)
In [51]:
y_true = df2[['y']]
y_pred = forecast1[['yhat']].iloc[:len(y_true)]
print(y_true.shape, y_pred.shape)
df_ypreds = pd.concat([y_true, y_pred],axis=1,sort=False).dropna()
df_ypreds.head().append(df_ypreds.tail())
Out[51]:
In [52]:
smape_np(df2['y'].to_numpy(), forecast1['yhat'].to_numpy())
Out[52]:
In [53]:
forecasts = [ 'forecast{}'.format(i) for i in range(1,7)]
smapes = [ smape(df2.y,eval('forecast{}.yhat'.format(i))) for i in range(1,7)]
smapes_fast = [ smape_fast(df2.y,eval('forecast{}.yhat'.format(i)))
for i in range(1,7)]
df_eval = pd.DataFrame({'forecast': forecasts,
'smape': smapes,
'fast_smape': smapes_fast})
df_eval
Out[53]:
Prophet Cross validations¶
In [54]:
from fbprophet.diagnostics import cross_validation
In [55]:
cv_results = cross_validation(m6,
initial='360 days',
period='30 days',
horizon='60 days')
In [56]:
smape_cv = smape(cv_results.y, cv_results.yhat)
smape_cv
Out[56]: